#include <cstdio>
#include <cstdlib>
#include <tuple>
#include <cctype>
#include <type_traits>
#include "cell.h"
#include "elemtag.hpp"
#include "esmd_types.h"
#include "hmini.hpp"
#include "list_cpp.hpp"
#include "htab_cpp.hpp"
#include "listed.hpp"
#include "memory_cpp.hpp"
#include "param_mapper.hpp"
#include "tuple_hash.hpp"
#include "unitconv.hpp"
#include "utils.hpp"
#include "scanner.hpp"
#include "treduce.hpp"
#include "gmxtop.hpp"
#include "scanner.hpp"
#include "listed_impl.hpp"

const int MAX_LINE_LEN = 65536;
using gmxscanner = scanner<MAX_LINE_LEN>;

template <typename... Us>
bool scan_param(gmxvparam<Us...> &vp, gmxscanner &scan, int ifunc, gmxfunclist<> flist) {
  printf("Unsupported function code %d\n", ifunc);
  return false;
}
template <typename T, typename... Ts, typename... Us>
bool scan_param(gmxvparam<Us...> &vp, gmxscanner &scan, int ifunc, gmxfunclist<T, Ts...> flist) {
  // puts(__PRETTY_FUNCTION__);
  if (ifunc == T::ifunc) {
    return scan.scan_all(vp.template as<typename T::param_t>().tie());
  } else {
    return scan_param(vp, scan, ifunc, gmxfunclist<Ts...>{});
  }
}
template <typename... Us>
bool scan_param(gmxvparam<Us...> &vp, gmxscanner &scan, int ifunc) {
  return scan_param(vp, scan, ifunc, gmxfunclist<Us...>{});
}
double convert_unit_single(double value, unit from, unit lunit, unit eunit, unit runit) {
  //e = kg.m.s-2
  unit r0 = from / eunit(from.p_mass); //Firstly use energy term to elimitate mass terms
  unit r1 = r0 / lunit(r0.p_leng);     //then eliminate length and rad terms
  unit r2 = r1 / runit(r1.p_rad);
  return value * r2.value;
}
template <typename... Us>
void convert_unit(gmxvparam<Us...> &vp, int ifunc, unit lunit, unit eunit, unit runit, gmxfunclist<> flist) {
  printf("Unsupported function code %d\n", ifunc);
}
template <typename... T1s, typename... T2s, int... Is>
auto convert_unit_seq(std::tuple<T1s...> values, std::tuple<T2s...> units, utils::iseq<int, Is...> seq, unit lunit, unit eunit, unit runit) {
  return std::make_tuple(convert_unit_single(std::get<Is>(values), std::get<Is>(units), lunit, eunit, runit)...);
}
template <typename T, typename... Ts, typename... Us>
void convert_unit(gmxvparam<Us...> &vp, int ifunc, unit lunit, unit eunit, unit runit, gmxfunclist<T, Ts...> flist) {
  if (ifunc == T::ifunc) {
    auto params = vp.template as<typename T::param_t>().tie();
    auto units = T::param_t::units_gmx;
    params = convert_unit_seq(params, units, utils::make_iseq<int, std::tuple_size<decltype(units)>{}>{}, lunit, eunit, runit);
  } else {
    convert_unit(vp, ifunc, lunit, eunit, runit, gmxfunclist<Ts...>{});
  }
}
template <typename... Us>
void convert_unit(gmxvparam<Us...> &vp, int ifunc, unit lunit, unit eunit, unit runit) {
  convert_unit(vp, ifunc, lunit, eunit, runit, gmxfunclist<Us...>{});
}
template <typename T>
bool scan_top(T &top, gmxscanner &sc) {
  if (!sc.scan_all(top.tie_tags()))
    return false;
  top.ifunc = 0;

  if (sc.scan_all(top.ifunc)) {
    top.override_param = scan_param(top.param, sc, top.ifunc);
  }
  //if (std::is_same<T, gmxtop::pair>::value) printf("ifunc: %d\n", top.ifunc);
  return true;
}
template <typename T>
bool scan_type(T &type, gmxscanner &sc) {
  return sc.scan_all(type.tie()) && scan_param(type.param, sc, type.ifunc);
}
static char *fgets_gmxtop(char *line, FILE *f) {
  char *p = line;
  if (!fgets_unlocked(p, MAX_LINE_LEN, f))
    return NULL;
  char last = 0;
  while ((*p != '\n') & (*p != 0))
    last = *(p++);
  while (last == '\\') {
    p--;
    if (!fgets_unlocked(p, MAX_LINE_LEN, f))
      break;
    while ((*p != '\n') & (*p != 0))
      last = *(p++);
  }
  p = line;
  while ((*p != ';') & (*p != 0))
    p++;
  *p = 0;
  return line;
}
gmxtop::gmxtop() : bond_types("bond_types"),
                   angle_types("angle_types"),
                   dihedral_types("dihedral_types"),
                   dirtab("directive htab"),
                   pair_types("pair_types"),
                   constraint_types("constraint types"),
                   molecule_types("molecule_types"),
                   molecules("molecules"),
                   atom_types("atom_types")
//  bond_mapper("bond_mapper"),
//  angle_mapper("angle_mapper"),
//  dihedral_mapper("dihedral_mapper")
{
  dirtab["defaults"] = DIR_DEFAULTS;
  dirtab["atomtypes"] = DIR_ATOMTYPES;
  dirtab["bondtypes"] = DIR_BONDTYPES;
  dirtab["pairtypes"] = DIR_PAIRTYPES;
  dirtab["angletypes"] = DIR_ANGLETYPES;
  dirtab["dihedraltypes"] = DIR_DIHEDRALTYPES;
  dirtab["constrainttypes"] = DIR_CONSTRAINTTYPES;
  dirtab["nonbonded_params"] = DIR_NONBONDED_PARAMS;
  dirtab["moleculetype"] = DIR_MOLECULETYPE;
  dirtab["atoms"] = DIR_ATOMS;
  dirtab["bonds"] = DIR_BONDS;
  dirtab["pairs"] = DIR_PAIRS;
  dirtab["angles"] = DIR_ANGLES;
  dirtab["constraints"] = DIR_CONSTRAINTS;
  dirtab["dihedrals"] = DIR_DIHEDRALS;
  dirtab["exclusions"] = DIR_EXCLUSIONS;
  dirtab["system"] = DIR_SYSTEM;
  dirtab["molecules"] = DIR_MOLECULES;
  dirtab["settles"] = DIR_SETTLES;
}
gmxtop::directive gmxtop::parse_directive(const char *line) {
  const char *p = line;
  while ((*p != 0) & (*p != '['))
    p++;
  if (*p == 0)
    return DIR_BROKEN;
  p++;
  while ((*p != 0) & (isspace(*p) != 0))
    p++;
  if (*p == 0)
    return DIR_BROKEN;
  const char *s = p;
  while ((*p != 0) & (*p != ']'))
    p++;
  if (*p == 0)
    return DIR_BROKEN;
  while (isspace(*(p - 1))) {
    p--;
  }
  char directive[256];
  strncpy(directive, s, p - s);
  directive[p - s] = 0;

  if (dirtab.contains(directive))
    return dirtab[directive];
  else {
    if (verbose) printf("directive [ %s ] not classified\n", directive);
    return DIR_OTHERS;
  }
}
gmxtop::directive gmxtop::parse_others(FILE *f) {
  char line[MAX_LINE_LEN];
  gmxscanner tz;
  while (fgets_gmxtop(line, f)) {
    if (tz.set_input(line) == '[')
      return parse_directive(line);
  }
  return DIR_FILEEND;
}

gmxtop::directive gmxtop::parse_moleculetype(FILE *f) {
  char line[MAX_LINE_LEN];
  gmxscanner tz;
  int cnt = 0;
  molecule_type *mt = new (molecule_types.get_slot()) molecule_type();
  while (fgets_gmxtop(line, f)) {
    char first = tz.set_input(line);
    if (first == '[')
      return parse_directive(line);

    if (tz.scan_all(mt->name, mt->nrexcl)) {
      if (cnt > 0) {
        if (verbose) printf("Warning: duplicate entry in moleculetype");
      }
      cnt++;
    }
  }
  return DIR_FILEEND;
}
template <typename T>
gmxtop::directive gmxtop::parse_types(esmd::list<T> &list, FILE *f) {
  char line[MAX_LINE_LEN];
  gmxscanner tz;

  while (fgets_gmxtop(line, f)) {
    char first = tz.set_input(line);
    if (first == '[') {
      return parse_directive(line);
    }
    T param;
    if (scan_type(param, tz)) {
      list.append(param);
    }
  }
  return DIR_FILEEND;
}
template <typename T>
gmxtop::directive gmxtop::parse_topologies(esmd::list<T> &list, FILE *f) {
  char line[MAX_LINE_LEN];
  gmxscanner tz;
  while (fgets_gmxtop(line, f)) {
    char first = tz.set_input(line);
    if (first == '[') {
      return parse_directive(line);
    }
    T top;
    if (scan_top(top, tz)) {
      list.append(top);
    }
  }
  return DIR_FILEEND;
}
gmxtop::directive gmxtop::parse_system(FILE *f) {
  char line[MAX_LINE_LEN];
  gmxscanner tz;
  while (fgets_gmxtop(line, f)) {
    char first = tz.set_input(line);
    if (first == '[')
      return parse_directive(line);
    if (first != 0)
      systemname = strdup(line);
  }
  return DIR_FILEEND;
}

template <typename T>
gmxtop::directive gmxtop::parse_to_list(esmd::list<T> &list, FILE *f) {
  char line[MAX_LINE_LEN];
  gmxscanner tz;
  while (fgets_gmxtop(line, f)) {
    char first = tz.set_input(line);
    if (first == '[')
      return parse_directive(line);
    T dat;
    if (tz.scan_all(dat.tie()))
      list.append(dat);
  }
  return DIR_FILEEND;
}
gmxtop::directive gmxtop::parse_defaults(FILE *f) {
  char line[MAX_LINE_LEN];
  gmxscanner tz;
  while (fgets_gmxtop(line, f)) {
    char first = tz.set_input(line);
    if (first == '[')
      return parse_directive(line);
    tz.scan_all(defaults.tie());
  }
  return DIR_FILEEND;
}

gmxtop::directive gmxtop::parse_exclusions(esmd::list<gmxtop::exclusion> &list, FILE *f) {
  char line[MAX_LINE_LEN];
  gmxscanner tz;
  while (fgets_gmxtop(line, f)) {
    char first = tz.set_input(line);
    if (first == '[')
      return parse_directive(line);
    exclusion e;
    if (!tz.scan_all(e.i))
      continue;
    while (tz.scan_all(e.j))
      list.append(e);
  }
  return DIR_FILEEND;
}

gmxtop::directive gmxtop::parse_settles(gmxtop::settle &set, FILE *f){
  char line[MAX_LINE_LEN];
  gmxscanner tz;
  while (fgets_gmxtop(line, f)) {
    char first = tz.set_input(line);
    if (first == '[')
      return parse_directive(line);
    if (!set.present)
      set.present = tz.scan_all(set.i, set.j, set.r0oh, set.r0hh);
  }
  return DIR_FILEEND;
}

void gmxtop::parse_file(FILE *f) {
  directive next_section = DIR_OTHERS;
  while (next_section != DIR_FILEEND && next_section != DIR_BROKEN) {
    switch (next_section) {
    case DIR_DEFAULTS:
      next_section = parse_defaults(f);
      break;
    case DIR_ATOMTYPES:
      next_section = parse_to_list(atom_types, f);
      break;
    case DIR_BONDTYPES:
      next_section = parse_types(bond_types, f);
      break;
    case DIR_ANGLETYPES:
      next_section = parse_types(angle_types, f);
      break;
    case DIR_DIHEDRALTYPES:
      next_section = parse_types(dihedral_types, f);
      break;
    case DIR_PAIRTYPES:
      next_section = parse_types(pair_types, f);
      break;
    case DIR_CONSTRAINTTYPES:
      next_section = parse_types(constraint_types, f);
      break;
    case DIR_MOLECULETYPE:
      next_section = parse_moleculetype(f);
      break;
    case DIR_ATOMS:
      next_section = molecule_types.size() ? parse_to_list(molecule_types[molecule_types.size() - 1].atoms, f) : parse_others(f);
      break;
    case DIR_BONDS:
      next_section = molecule_types.size() ? parse_topologies(molecule_types[molecule_types.size() - 1].bonds, f) : parse_others(f);
      break;
    case DIR_ANGLES:
      next_section = molecule_types.size() ? parse_topologies(molecule_types[molecule_types.size() - 1].angles, f) : parse_others(f);
      break;
    case DIR_DIHEDRALS:
      next_section = molecule_types.size() ? parse_topologies(molecule_types[molecule_types.size() - 1].dihedrals, f) : parse_others(f);
      break;
    case DIR_CONSTRAINTS:
      next_section = molecule_types.size() ? parse_topologies(molecule_types[molecule_types.size() - 1].constraints, f) : parse_others(f);
      break;
    case DIR_PAIRS:
      next_section = molecule_types.size() ? parse_topologies(molecule_types[molecule_types.size() - 1].pairs, f) : parse_others(f);
      break;
    case DIR_EXCLUSIONS:
      next_section = molecule_types.size() ? parse_exclusions(molecule_types[molecule_types.size() - 1].excls, f) : parse_others(f);
      break;
    case DIR_SETTLES:
      next_section = molecule_types.size() ? parse_settles(molecule_types[molecule_types.size() - 1].settles, f) : parse_others(f);
      break;
    case DIR_SYSTEM:
      next_section = parse_system(f);
      break;
    case DIR_MOLECULES:
      next_section = parse_to_list(molecules, f);
      break;
    default:
      if (verbose) printf("directive code %d not classified\n", next_section);
      next_section = parse_others(f);
    }
  }
  bond_types.sort();
  angle_types.sort();
  dihedral_types.sort();
  pair_types.sort();
  if (verbose) printf("%d\n", next_section);
}

template <typename T>
int print_t(T t) {
  if (std::is_same<T, int>::value) {
    printf("%d ", t);
  } else if (std::is_same<T, double>::value) {
    printf("%f ", t);
  } else if (std::is_same<T, char *>::value) {
    printf("%s ", t);
  } else {
    printf("UNSUPPORTED_TYPE ");
    return 0;
  }
  return 1;
}
int print_t(elemtag tag) {
  return printf("%s ", tag.str) == 1;
}
template <typename... Ts, int... Is>
void print_tuple(std::tuple<Ts...> tpl, utils::iseq<int, Is...> seq) {
  int a[] = {print_t(std::get<Is>(tpl))...};
}
template <typename... Ts>
void print_tuple(std::tuple<Ts...> tpl) {
  print_tuple(tpl, utils::make_iseq<int, sizeof...(Ts)>{});
}
template <typename... Us>
void print_param(gmxvparam<Us...> vp, int ifunc, gmxfunclist<> flist) {
  printf("UNSUPPORTED_IFUNC%d\n", ifunc);
}
template <typename... Us, typename T, typename... Ts>
void print_param(gmxvparam<Us...> vp, int ifunc, gmxfunclist<T, Ts...> flist) {
  // puts(__PRETTY_FUNCTION__);
  if (T::ifunc == ifunc) {
    print_tuple(vp.template as<typename T::param_t>().tie());
  } else {
    print_param(vp, ifunc, gmxfunclist<Ts...>{});
  }
}
template <typename... Us>
void print_param(gmxvparam<Us...> vp, int ifunc) {
  print_param(vp, ifunc, gmxfunclist<Us...>{});
}
void gmxtop::unitconv(unit lunit, unit eunit, unit runit) {
  for (int i = 0; i < bond_types.size(); i++) {
    bond_type &p = bond_types[i];
    convert_unit(p.param, p.ifunc, lunit, eunit, runit);
  }
  for (int i = 0; i < angle_types.size(); i++) {
    angle_type &p = angle_types[i];
    convert_unit(p.param, p.ifunc, lunit, eunit, runit);
  }
  for (int i = 0; i < dihedral_types.size(); i++) {
    dihedral_type &p = dihedral_types[i];
    convert_unit(p.param, p.ifunc, lunit, eunit, runit);
  }
  for (auto &moltype : molecule_types){
    for (auto &b : moltype.bonds){
      if (b.override_param){
        convert_unit(b.param, b.ifunc, lunit, eunit, runit);
      }
    }
    for (auto &a : moltype.angles){
      if (a.override_param){
        convert_unit(a.param, a.ifunc, lunit, eunit, runit);
      }
    }
    for (auto &d : moltype.dihedrals){
      if (d.override_param){
        convert_unit(d.param, d.ifunc, lunit, eunit, runit);
      }
    }
  }
  // auto factor_sig = (pow(2, 1.0/6.0) * units::nm / lunit).value;
  // auto factor_eps = (units::kJ / units::mole / eunit).value;
  // printf("%g %g\n", factor_sig, factor_eps);
  // for (auto &atype : atom_types) {
  //   atype.sigma *= factor_sig;
  //   atype.epsilon *= factor_eps;
  // }
}
template <typename T>
void print_type(T t) {
  print_tuple(t.tie());
  print_param(t.param, t.ifunc);
  puts("");
}

template <typename T, typename P>
void find_parameters_list(esmd::list<T> &top_list, const gmxtop::molecule_type &mol, P &mapper) {
  for (int i = 0; i < top_list.size(); i++) {
    T &top = top_list[i];
    if (top.ifunc != 0) {
      if (top.override_param) {
        top.pid = mapper.map(top.ifunc, top.param);
      } else {
        for (int k = 0; k < T::nord; k++) {
          top.pid = mapper.get(mol.get_types(top.tpl_tags(k)), top.ifunc);
          if (top.pid != -1)
            break;
        }
      }
    } else {
      for (int k = 0; k < T::nord; k++) {
        std::tie(top.ifunc, top.pid) = mapper.get(mol.get_types(top.tpl_tags(k)));
        if (top.pid != -1)
          break;
      }
    }
  }
}

void summary(gmxtop &top) {
  if (top.verbose) printf("topology has %ld atom types, %ld bond types, %ld angle types, %ld dihedral types, %ld pair types\n", top.atom_types.size(), top.bond_types.size(), top.angle_types.size(), top.dihedral_types.size(), top.pair_types.size());
  if (top.verbose) printf("topology has %ld molecule types", top.molecule_types.size());
  for (int i = 0; i < top.molecule_types.size(); i++) {
    gmxtop::molecule_type *mt = &top.molecule_types[i];
    if (top.verbose) printf("molecule %d's name is %s\n", i, mt->name);
    if (top.verbose) printf("molecule %d has %ld atoms, %ld bonds, %ld angles, %ld dihedrals\n", i, mt->atoms.size(), mt->bonds.size(), mt->angles.size(), mt->dihedrals.size());
  }
  if (top.verbose) printf("system name is: %s\n", top.systemname);
  for (int i = 0; i < top.molecules.size(); i++) {
    gmxtop::molecule &mol = top.molecules[i];
    if (top.verbose) printf("molecule: %s count: %ld\n", mol.name, mol.cnt);
  }
}

__always_inline bool bit_in_set(int bit, size_t set){
  return (1L << bit) & set;
}
__always_inline size_t make_bitset(){
  return 0;
}
template<typename Arg, typename ...Args>
__always_inline size_t make_bitset(Arg first, Args... others){
  return (1L << first) | make_bitset(others...);
}
template<typename TopoType, typename ParamType>
struct gmx_adapter{
  size_t begin, end, i;
  int ifunc_sels;
  esmd::list<TopoType> &param_list;
  ParamType param;
  gmx_adapter(esmd::list<TopoType> &param_list, int ifunc_sels) : param_list(param_list), ifunc_sels(ifunc_sels) {}
  void set(size_t begin, size_t end) {
    this->begin = begin;
    this->end = end;
    this->i = -1;
  }
  void reset() {
    i = -1;
  }
  bool shift(){
    while (begin < end) {
      if (bit_in_set(param_list[begin].ifunc, ifunc_sels)) {
        i = begin;
        begin ++;
        return true;
      }
      begin ++;
    }
    return false;
  }
  bool has_data(){
    return i != -1;
  }
};
struct gmx_angle_to_ub : gmx_adapter<gmxtop::angle_type, urey_bradly_param> {
  gmx_angle_to_ub(esmd::list<gmxtop::angle_type> &param_list, int ifunc_sels) : gmx_adapter(param_list, ifunc_sels) {}
  bool next(){
    if (shift()){
      auto &param_angle = param_list[i].param.template as<harmonic_angle_param>();
      param = {param_angle.theta0, param_angle.ktheta, 0, 0};
      return true;
    }
    return false;
  }
};

struct gmx_pair_route : gmx_adapter<gmxtop::pair_type, routed_ljcoul_param> {
  gmx_pair_route(esmd::list<gmxtop::pair_type> &param_list, int ifunc_sels) : gmx_adapter(param_list, ifunc_sels){}
  bool next(){
    if (shift()){
      //printf("%d\n", i);
      auto &param_lj = param_list[i].param.template as<ljcoul_param>();
      this->param = {param_lj.sigma, param_lj.epsilon};
      return true;
    }
    return false;
  }
};

template<typename FromParam, typename ToParam>
struct gmx_pass : public gmx_adapter<FromParam, ToParam> {
  gmx_pass(esmd::list<FromParam> &param_list, int ifunc_sels) : gmx_adapter<FromParam, ToParam>(param_list, ifunc_sels) {}
  bool next(){
    if (this->shift()){
      this->param.tie() = this->param_list[this->i].param.template as<ToParam>().tie();
      return true;
    }
    return false;
  }
};
struct gmx_cosine_to_rb : public gmx_adapter<gmxtop::dihedral_type, extended_ryckaert_bellmans_param> {
  gmx_cosine_to_rb(esmd::list<gmxtop::dihedral_type> &param_list, int ifunc_sels) : gmx_adapter(param_list, ifunc_sels) {}

  bool next(){
    reset();
    param.tie() = std::make_tuple(0, 0, 0, 0, 0, 0, 0);
    while (shift()) {
      if (!bit_in_set(param_list[i].ifunc, ifunc_sels)) continue;
      auto param_cos = (cosine_torsion_param)param_list[i].param;
      param += param_cos;
    }
    return has_data();
  }
};

template<typename ParamAdapter, typename ToParam>
void gmxtop::export_listed_adapt(listed_force<ToParam> &listed, param_mapper<ToParam> &mapper, int ifunc_sels){
  typedef gmx_top_of<ToParam> TopoType;
  tagint atom_off = -1;
  auto &top_types = this->get_params<typename TopoType::param_type>();

  ParamAdapter adapt(top_types, ifunc_sels);
  for (int imol = 0; imol < molecules.size(); imol ++){
    molecule &mol = molecules[imol];
    molecule_type &moltype = get_moltype(mol.name);
    size_t list_start = listed.topo.size();
    tagint natom = moltype.atoms.size();
    auto &tops = moltype.get_tops<TopoType>();
    //printf("%d\n", tops.size());
    typedef typename TopoType::param_type ParamType;
    //printf("%ld\n", atom_off);
    for (int itop = 0; itop < tops.size(); itop ++){
      auto &topo = tops[itop];
      if (topo.ifunc != 0 && !bit_in_set(topo.ifunc, ifunc_sels)) continue;
      auto tags = topo.tie_tags();

      if (topo.override_param){
        auto entry = listed_force_tag_entry<ToParam>::from_tuple(tags, mapper.map(topo.param.template as<ToParam>()));
        entry.shift(atom_off);
        listed.topo.append(entry);
      } else {
        for (int ord = 0; ord < TopoType::nord; ord ++){
          auto types = moltype.get_types(topo.tpl_tags(ord));
          int ifunc_begin = topo.ifunc ? topo.ifunc : INT32_MIN;
          int ifunc_end = topo.ifunc ? topo.ifunc : INT32_MAX;
          ParamType key_begin, key_end;
          key_begin.tie() = std::tuple_cat(types, std::make_tuple(ifunc_begin));
          key_end.tie() = std::tuple_cat(types, std::make_tuple(ifunc_end));
          size_t begin = top_types.bisect_begin(key_begin);
          size_t end = top_types.bisect_end(key_end);
          if (begin < end) {
            adapt.set(begin, end);
            while (adapt.next()) {
              auto entry = listed_force_tag_entry<ToParam>::from_tuple(tags, mapper.map(adapt.param));
              entry.shift(atom_off);
              listed.topo.append(entry);
            }
            if (adapt.has_data()) break;
          }
        }
      }
    }
    size_t list_end = listed.topo.size();
    for (int icnt = 1; icnt < mol.cnt; icnt ++){
      for (int ilist = list_start; ilist < list_end; ilist ++) {
        auto entry = listed.topo[ilist];
        entry.shift(natom * icnt);
        listed.topo.append(entry);
      }
    }
    atom_off += mol.cnt * natom;
  }

  listed.param = mapper.param_list.extract();
  listed.nparam = mapper.param_list.size();
}

template<typename ToParam>
void gmxtop::export_listed(listed_force<ToParam> &listed, param_mapper<ToParam> &mapper, int ifunc_sels){
  typedef gmx_top_of<ToParam> TopoType;
  this->export_listed_adapt<gmx_pass<typename TopoType::param_type, ToParam>, ToParam>(listed, mapper, ifunc_sels);
}
template void gmxtop::export_listed(listed_force<harmonic_bond_param> &, param_mapper<harmonic_bond_param> &, int);
template void gmxtop::export_listed(listed_force<harmonic_angle_param> &, param_mapper<harmonic_angle_param> &, int);
template void gmxtop::export_listed(listed_force<urey_bradly_param> &, param_mapper<urey_bradly_param> &, int);
template void gmxtop::export_listed(listed_force<cosine_torsion_param> &, param_mapper<cosine_torsion_param> &, int);
template void gmxtop::export_listed(listed_force<ryckaert_bellmans_param> &, param_mapper<ryckaert_bellmans_param> &, int);
template void gmxtop::export_listed(listed_force<harmonic_improper_param> &, param_mapper<harmonic_improper_param> &,int);
template void gmxtop::export_listed_adapt<gmx_pair_route>(listed_force<routed_ljcoul_param> &, param_mapper<routed_ljcoul_param> &, int);

void gmxtop::build_exclusions(){
  for (int i = 0; i < molecule_types.size(); i ++){
    molecule_types[i].build_exclusions(*this);
  }
}
void gmxtop::molecule_type::build_exclusions(gmxtop &top){
  struct atom_pair{
    tagint src, dst;
    bool operator<(const atom_pair &o){
      return src < o.src;
    }
  };
  esmd::list<atom_pair> edges("build exclusions/temp bonds");
  for (tagint i = 0; i < bonds.size(); i ++){
    edges.append(atom_pair{bonds[i].i, bonds[i].j});
    edges.append(atom_pair{bonds[i].j, bonds[i].i});
  }
  esmd::htab<tuple_key<elemtag, elemtag>, bool> constraint_bonded("build exclusions/constraint ifunc");
  for (int i = 0; i < top.constraint_types.size(); i ++){
    if (top.constraint_types[i].ifunc == 1){
      constraint_bonded[std::make_tuple(top.constraint_types[i].i, top.constraint_types[i].j)] = true;
      constraint_bonded[std::make_tuple(top.constraint_types[i].j, top.constraint_types[i].i)] = true;
    }
  }
  for (tagint i = 0; i < constraints.size(); i ++){
    elemtag itype = atoms[constraints[i].i - 1].atom_name;
    elemtag jtype = atoms[constraints[i].j - 1].atom_name;
    if (constraints[i].ifunc == 1 || constraint_bonded.contains(std::make_tuple(itype, jtype))){
      edges.append(atom_pair{constraints[i].i, constraints[i].j});
      edges.append(atom_pair{constraints[i].j, constraints[i].i});
    }
  }

  edges.sort();
  tagint *first_edge = esmd::allocate<tagint>(atoms.size() + 2, "build exclusions/temp heads");
  tagint atom_ptr = 0;
  for (tagint i = 0; i < edges.size(); i ++){
    while (atom_ptr < edges[i].src) {
      atom_ptr ++;
      first_edge[atom_ptr] = i;
    }
  }
  while (atom_ptr <= atoms.size()){
    atom_ptr ++;
    first_edge[atom_ptr] = edges.size();
  }
  esmd::list<tagint> queue("build exclusions/queue");
  for (tagint i = 1; i <= atoms.size(); i ++){
    mini_hset<tagint, 127> excl_set;
    excl_set.insert(i);
    queue.clear();
    queue.append(i);
    int head = 0;
    for (int dist = 0; dist < nrexcl; dist ++){
      int tail = queue.size();
      for (int jj = head; jj < tail; jj ++){
        tagint j = queue[jj];
        for (tagint kk = first_edge[j]; kk < first_edge[j + 1]; kk ++){
          int k = edges[kk].dst;
          if (!excl_set.contains(k)) {
            queue.append(k);
            excl_set.insert(k);
          }
        }
      }
      head = tail;
    }
    for (int jj = 1; jj < queue.size(); jj ++){
      excls.append(exclusion{i, queue[jj]});
    }
  }
  // if (!strcmp(name, "water1")){
  //   printf("exclusions in water1:");
  //   for (auto excl: excls){
  //     printf("excl %ld %ld\n", excl.i, excl.j);
  //   }
  // }
  esmd::deallocate(first_edge);
}

void gmxtop::molecule_type::sort_rigids() {
  for (auto &b : bonds){
    if (atoms[b.i - 1].mass < atoms[b.j - 1].mass) {
      std::swap(b.i, b.j);
    }
  }
}
void gmxtop::sort_rigids() {
  for (auto &moltype : molecule_types){
    moltype.sort_rigids();
  }
}
void gmxtop::molecule_type::find_clusters(listed_force<rigid_param> &rigids, param_mapper<rigid_param> &param_mapper, esmd::htab<tuple_key<elemtag, elemtag>, double> &r0s) {
  auto is_hydrogen = [](atom &a){
    return (a.type.str[0] == 'H' || a.type.str[0] == 'h');
  };
  int cnt = 0;
  for (auto &b : bonds){
    atom &iatom = atoms[b.i - 1];
    atom &jatom = atoms[b.j - 1];
    if (is_hydrogen(iatom) || is_hydrogen(jatom)) {
      cnt ++;
      constraint c{b.j, b.i, 1};
      if (b.override_param){
        c.override_param = true;
        c.param.as<distance_constraint>().r0 = b.param.as<harmonic_bond_param>().r0;
      }
      constraints.append(c);
    }
  }
  int bond_cnt = 0;
  /*filter out hydrogen related bonds*/
  for (auto b : bonds){
    atom &iatom = atoms[b.i - 1];
    atom &jatom = atoms[b.j - 1];
    if (!is_hydrogen(iatom) && !is_hydrogen(jatom)) {
      bonds[bond_cnt] = b;
      bond_cnt ++;
    }
  }
  bonds.set_size(bond_cnt);
  //printf("%d\n", cnt);
  for (auto &c : constraints){
    atom &iatom = atoms[c.i - 1];
    atom &jatom = atoms[c.j - 1];
    if (iatom.mass < jatom.mass) {
      tagint tmp = c.i;
      c.i = c.j;
      c.j = tmp;
    }
  }
  constraints.sort();
  //param_mapper<rigid_param> rmapper("rigid param mapper");
  //esmd::list<listed_force_entry<rigid_param, tagint>> rigids("rigid clusters");
  //if (verbose) printf("#constraints %ld\n", constraints.size());
  for (int i = 0; i < constraints.size(); ){
    listed_force_entry<rigid_param, tagint> rigid;

    rigid.id[0] = constraints[i].i;
    rigid_param par = {0, 0, 0, 0};
    int j = i;

    while (j < constraints.size() && constraints[j].i == constraints[i].i) {
      rigid.id[j - i + 1] = constraints[j].j;
      real r0;
      if (constraints[j].override_param) {
        r0 = constraints[j].param.as<distance_constraint>().r0;
      } else {
        assert(r0s.contains(std::make_tuple(atoms[constraints[j].i - 1].type, atoms[constraints[j].j - 1].type)));
        r0 = r0s[std::make_tuple(atoms[constraints[j].i - 1].type, atoms[constraints[j].j - 1].type)];
      }
      par.r0[j - i] = r0;
      j ++;
    }

    par.type = j - i + 1;

    for (int j = par.type; j < 4; j ++) rigid.id[j] = rigid.id[0];
    for (int jj = 1; jj < par.type; jj ++){
      for (int j = 1; j < par.type - 1; j ++){
        if (par.r0[j] < par.r0[j - 1]) {
          std::swap(par.r0[j], par.r0[j - 1]);
          std::swap(rigid.id[j], rigid.id[j + 1]);
        }
      }
    }
    rigid.pid = param_mapper.map(par);
    rigids.topo.append(rigid);
    //rigids.append(rigid);
    //puts("");
    i = j;
  }
  //rigids.nparam = param_mapper.param_list.size();
  // for (auto &p : rmapper.param_list){
  //   for (int j = 0; j < p.type - 1; j ++){
  //     printf("%f ", p.r0[j]);
  //   }
  //   puts("");
  // }
  // printf("#clusters %ld, #cluster types %ld\n", rigids.topo.size(), rmapper.param_list.size());
}

void gmxtop::find_clusters(listed_force<rigid_param> &rigids, param_mapper<rigid_param> &param_mapper){
  esmd::htab<tuple_key<elemtag, elemtag>, real> dist_tab("rigid disttab");
  for (int i = 0; i < bond_types.size(); i ++){
    dist_tab[std::make_tuple(bond_types[i].i, bond_types[i].j)] = bond_types[i].param.as<harmonic_bond_param>().std_len();
    dist_tab[std::make_tuple(bond_types[i].j, bond_types[i].i)] = bond_types[i].param.as<harmonic_bond_param>().std_len();
  }
  for (int i = 0; i < constraint_types.size(); i ++){
    dist_tab[std::make_tuple(constraint_types[i].i, constraint_types[i].j)] = constraint_types[i].param.as<distance_constraint>().r0;
    dist_tab[std::make_tuple(constraint_types[i].j, constraint_types[i].i)] = constraint_types[i].param.as<distance_constraint>().r0;
  }
  // for (int i = 0; i < molecule_types.size(); i ++){
  //   molecule_types[i].find_clusters(rigids, param_mapper, dist_tab);
  // }
  tagint offset = -1;
  for (auto &mol : molecules){
    auto &moltype = get_moltype(mol.name);
    size_t cur = rigids.topo.size();
    moltype.find_clusters(rigids, param_mapper, dist_tab);
    size_t cnt = rigids.topo.size() - cur;
    for (size_t i = 1; i < mol.cnt; i ++){
      for (size_t j = 0; j < cnt; j ++){
        auto rig = rigids.topo[cur + j];
        rigids.topo.append(rig);
      }
    }
    for (size_t i = 0; i < mol.cnt; i ++){
      for (size_t j = 0; j < cnt; j ++){
        rigids.topo[cur + i * cnt + j].shift(offset);
      }
      offset += moltype.atoms.size();
    }
  }
  rigids.param = param_mapper.param_list.extract();
  rigids.nparam = param_mapper.param_list.size();
}

void gmxtop::molecule_type::route_pairs(){
  esmd::htab<tuple_key<tagint, tagint>, tagint> routes("route pairs/route map");
  for (auto &d : dihedrals){
    if ((1 << d.ifunc & dihedral::impr_ifuncs) == 0) {
      routes[std::make_tuple(d.i, d.l)] = d.j;
      routes[std::make_tuple(d.l, d.i)] = d.j;
    }
  }

  for (auto &p : pairs) {
    auto key = std::make_tuple(p.i, p.j);
    assert(routes.contains(key));
    routed_pair rp;
    rp.i = p.i;
    rp.j = p.j;
    rp.k = routes[key];
    //printf("ifunc: %d\n", p.ifunc);
    rp.ifunc = p.ifunc;
    rp.override_param = p.override_param;
    rp.param = p.param;
    routed_pairs.append(rp);
  }
  //if (verbose) printf("%d %d %d\n", pairs.size(), routed_pairs.size(), routed_pairs[0].ifunc);
}

void gmxtop::route_pairs(){
  for (auto &mol : molecule_types) {
    mol.route_pairs();
  }
}

void gmxtop::fudge_pairtypes() {
  esmd::htab<tuple_key<elemtag, elemtag>, ljcoul_param> params("param type");
  esmd::htab<tag_key, atom_type> atom_map("atom map");
  for (auto &pt : pair_types) {
    ljcoul_param param = pt.param.as<ljcoul_param>();
    params[std::make_tuple(pt.i, pt.j)] = param;
  }
  for (auto &it : atom_types) {
    //atom_map[at.name.ival] = at;
    for (auto &jt : atom_types) {
      auto key = std::make_tuple(it.name, jt.name);
      if (!params.contains(key)) {
        ljcoul_param param;
        param.epsilon = defaults.fudge_lj * sqrt(it.epsilon * jt.epsilon);
        param.sigma = 0.5 * (it.sigma * jt.sigma);
        param.qq = defaults.fugde_qq * it.charge * jt.charge;
        params[key] = param;
      }
    }
  }
  pair_types.clear();
  for (auto &mp : params) {
    pair_type pt;
    std::tie(pt.i, pt.j) = mp.key;
    pt.ifunc = 1;
    pt.param.as<ljcoul_param>() = mp.val;
    pair_types.append(pt);
  }
  pair_types.sort();
}
void gmxtop::create_vdw14_params(esmd::htab<tag_key, std::tuple<real, real, real, real>> &param_map){
  for (auto &at : atom_types){
    real eps = at.epsilon;
    real sig = at.sigma;
    real eps14 = at.epsilon * defaults.fudge_lj;
    real sig14 = at.sigma;
    for (auto &pt : pair_types){
      if (pt.i.ival == at.name.ival && pt.j.ival == at.name.ival) {
        eps14 = pt.param.as<ljcoul_param>().epsilon;
        sig14 = pt.param.as<ljcoul_param>().sigma;
      }
    }
    param_map[at.name.ival] = std::make_tuple(eps, sig, eps14, sig14);
  }
}
void gmxtop::export_excl_scal(cellgrid_t &grid) {
  esmd::list<exclusion> bidir_excls("gmxtop/bidir_excls");
  esmd::list<exclusion> bidir_scals("gmxtop/bidir_scals");
  tagint offset = -1;
  for (auto &mol : molecules) {
    auto &moltype = get_moltype(mol.name);
    for (int i = 0; i < mol.cnt; i ++){
      for (auto &excl : moltype.excls){
        bidir_excls.append(exclusion{excl.i + offset, excl.j + offset});
      }
      for (auto &pair : moltype.pairs){
        bidir_scals.append(exclusion{pair.i + offset, pair.j + offset});
        bidir_scals.append(exclusion{pair.j + offset, pair.i + offset});
      }
      offset += moltype.atoms.size();
    }
  }
  bidir_excls.sort();
  bidir_scals.sort();
  if (verbose) printf("bidirectional exclusions: %d, scales: %d\n", bidir_excls.size(), bidir_scals.size());
}

void gmxtop::export_special_pairs(listed_force<special_pair> &excls) {
  typedef listed_force_entry<special_pair, long> special_entry;
  tagint offset = -1;
  for (auto &mol : molecules) {
    auto &moltype = get_moltype(mol.name);
    for (int i = 0; i < mol.cnt; i ++){
      for (auto &excl : moltype.excls){
        excls.topo.append(special_entry{excl.i + offset, excl.j + offset, special_pair::EXCL});
      }
      for (auto &pair : moltype.pairs){
        excls.topo.append(special_entry{pair.i + offset, pair.j + offset, special_pair::SCAL});
        excls.topo.append(special_entry{pair.j + offset, pair.i + offset, special_pair::SCAL});
      }
      offset += moltype.atoms.size();
    }
  }
  
}

void gmxtop::molecule_type::convert_settle(){
  if (settles.present){
    // printf("converting settles for %s\n", name);
    // printf("original:\n");
    // for (auto bond: bonds){
    //   printf("bond %ld %ld %d\n", bond.i, bond.j, bond.ifunc);
    // }

    // for (auto excl: excls){
    //   printf("excl %ld %ld\n", excl.i, excl.j);
    // }

    bonds.clear();
    excls.clear();    
    for (int i = 0; i < atoms.size(); i ++){
      if (settles.i != i + 1)
        bonds.append(bond{.i = settles.i, .j = i + 1, .ifunc = 1});
    }
    // printf("converted to:\n");
    // for (auto bond: bonds){
    //   printf("bond %ld %ld %d\n", bond.i, bond.j, bond.ifunc);
    // }

    // for (auto excl: excls){
    //   printf("excl %ld %ld\n", excl.i, excl.j);
    // }

  }
}

void gmxtop::convert_settle(){
  for (auto &moltype : molecule_types){
    moltype.convert_settle();
  }
}
